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Abstract . The National Ocean Survey Is devel- 
oping an automated system to derive parameters of 
horizontal crustal motion from existing geodetic 
data by the process of least-squares estimation. 
The estimated parameters will describe crustal 
motion as a function of geographic position. The 
system will first be tested in the Imperial Val- 
ley region of southern California, using data 
from 8 individual field projects spanning four 
decades of time. 


Introduction 

Global models for tectonic activity hypoth- 
esize the existence of rigid plates rotating with 
constant velocity. In contrast, local crustal 
motion as observed by geodetic and geophysical 
instrumentation varies from nearly continuous 
creep to the stop and go process associated with 
large earthquakes . To better understand the 
transition from global to local phenomena, the 
National Geodetic Survey of the National Ocean 
Survey is performing several studies of the 
geodetic data of the past 100 years on a re- 
gional level. These studies are designed to 
establish the pattern of horizontal crustal 
motion in areas from 100 to 300 km in diameter. 
This paper describes one of the techniques being 
used and some preliminary results obtained from 
a pilot study of the Imperial Valley area in 
southern California. 

The participation of the National Geodetic 
Survey in crustal motion study is required by 
the forthcoming redefinition of the North Amer- 
ican Horizontal Datum. A new adjustment of the 
entire U.S. control network will accompany this 
redefinition, and new positions will be pub- 
lished in 1983 for all stations of the con- 
trol network. The following arguments are pre- 
sented to support the geodetic community's need 
for a better understanding of crustal motion. 

1. In any adjustment which Incorporates data 
from different epochs, especially in an 
area of crustal movement, the observations 
need to be reduced to a model of the earth 
which allows geodetic positions to vary 
with time. 

2., A model for crustal motion is needed for 
predicting the changes in position of pub- 
lished stations. To the extent feasible, 
parameters of motion could be published 
in 1983 along with station positions in a 
fashion similar to star catalogs. 

3. A better understanding of crustal motion 
will help to better define requirements 
for reobserving disrupted sections of the 
control network. 
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The Technique 

In most studies of crustal motion the usual 
technique is to compare different surveys of the 
same geodetic network, two epochs at a time. 

The technique to be discussed here will differ 
out of necessity. The geographic areas of study 
will generally be larger than the area covered 
by any one field project, and in most cases the 
various field projects will only partially over- 
lap one another since most of them were observed 
to establish geodetic control where it pre- 
viously did not exist, not for crustal motion 
study. The basic technique is to estimate param- 
eters describing crustal motion by a simulta- 
neous least-squares adjustment of all the per- 
tinent geodetic data. This is accomplished by 
introducing into the adjustment a mathematical 
model which describes station positions as a 
function of time. The following paragraph 
describes the model which was used in the Impe- 
rial Valley pilot test. 

In the model the region of study can consist 
of one or more subregions. Existing fault lines 
will usually provide the boundaries between 
subregions. The latitude ij)(. and the longitude 
Xf. of a geodetic station in the ith subregion 
at time t are given by the formulas. 

♦t • * 'l.l + 'l.S "-‘o)' 

*t ■ + 'l,2 

Here is a fixed time of reference, and 
(*)>t geodetic coordinates of the 

station at time t(j. Each f^^ j for l<j<4 is a 
function over the variables and X^. . In 

the first applications of the technique these 
functions will be of the form 
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Note that this models the motion as a contin- 
uous function of time and a discontinuous func- 
tion of position with the discontinuities occur- 
ring along the boundaries between subregions. 
Existing horizontal survey data in the form of 
directions, distances, and azimuths can be input 
into the adjustment process to obtain the least- 
squares estimates for the unknown coordinates 
(<^t( 3 >^to) and the unknown coefficients 

This technique has several advantages over the 
standard technique of directly comparing two sets 
of measurements of the same quantities. It 
allows for the linking together of neighboring 
field projects into a single data set even though 
several years might exist between the times when 
the individual field projects were observed. It 
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allows for the rigorous Inclusion of astronomic 
azimuths which may have been observed separately 
and the inclusion of data at stations which have 
been destroyed. Additionally, the model pro- 
vides a built-in mechanism for interpolating the 
values of velocity and acceleration over the 
entire region of study. Finally, the model gen- 
eralizes the information contained in the data. 
This last point can be considered a disadvantage 
as well as an advantage. It is an advantage in 
so much as the concern is toward regional trends 
as opposed to local details. For example, local 
movement phenomena like hillside creep will be 
smoothed-out . On the other hand , it is a dis- 
advantage in that unmodeled variations in re- 
gional motion will also be smoothed-out. For 
example, motion across an unmodeled, yet active, 
fault will be interpreted as a continuous func- 
tion of position. For this reason the standard 
technique of directly comparing two sets of 
measurements over the same quantities will 
never be fully abandoned. Instead it will pro- 
vide the standards by which to evaluate the ac- 
curacy of the mathematical model. 

Although information is lost in the process 
of modeling the motion, a model is desired to 
provide a clear overall picture. When the in- 
adequacies of a model are identified, the model 
can be refined. The study of the Imperial 
Valley data was conducted as a pilot test to 
evaluate the above model, identify its inade- 
quacies, and suggest how the model might be 
changed to more accurately reflect the physi- 
cal situation. 



The Data 

Of the seismic areas in the United States, 
the Imperial Valley and its immediate sur- 
roundings represent the most frequently ob- . 
served part of the national network over any 
geographic area of comparable size. The basic 
network was observed in 1934. After the El 
Centro earthquake of May, 1940 (Richter scale 
magnitude =7.1) the network was repbserved 
in 1941 to determine a new set of positions for 
several geodetic stations. Observations of the 
basic network were performed in 1954-55 and 
again in 1967 for the specific purpose of study- 
ing post-seismic activity in the area. In the 
period 1974-76 two field projects along the 
southern extent of the network were observed as 
part of the Transcontinental Traverse. . This 
study also Includes two minor field projects 
in the area, a 1942 triangulation survey around 
the southern half of the Salton Sea and a 1959 
highway traverse survey extending parallel and 
about 10 km north of the Mexlco-Callfornia 
border. Figure 1 shows the essential part of 
the network which was studied and its relation- 
ship to the fault system in the area. 

Different parts of this data were previously 
investigated. The results of these investiga- 
tions can be used to evaluate the performance 
of the model. Displacement vectors for the 
region were reported by Meade [1948] for the 
1934-1941 period, Whitten [1956] for the 1941- 
55 period, and Gergen [1978] for the 1941-1976 
period. Miller et al. [1970] published dis- 






placement vectors and strain components for 
various periods using appropriate subsets of 
the 1934-41-55-67 data. A geophysical interpre- 
tation of the data was given by Scholz and Fitch 
[1969] by comparing the 1941-55 data to the dis- 
location model of Chinnery [1961]. Interpreta- 
tions of the data were also performed by Savage 
and Burford [1970], Barker [1976], and Thatcher 
[1978] . Each of these last three papers ana- 
lyzed the strain components derived from the 
1941-55-67 data by the method of Frank [1966]. 

Experiments with the Model 

The model as programmed for this pilot test 
allows the user to partition the region of study 
into three subregions and to solve for 15 coef- ‘ 
ficients (all terms in <(i and X up to degree four) 
for each of the four polynomials associated with 
a subregion. This gives the user a total of 180 
parameters with which to describe the motion. 

The experiments are to determine the appropri- 
ateness of these parameters. However, it was 
first necessary to establish the location of 
active faults. This was accomplished with geo- 
logical maps and an adjustment of the data to a 
model which does not Include any time param- 
eters. The ^geological maps located the known 
faults. The residuals obtained from the ad- 
justment identified which of these faults were 
the most active. In some Instances the maps 
were ambiguous as to the location of a station 
relative to the fault. Station B In figure 2 
illustrates the problem. These ambiguities were 
resolved by assuming one sense for the relative 
motion between opposite sides of the fault and 
checking whether the. angle o at B measured clock- 
wise from A to C Is increasing or decreasing with 
time. In figure 2 right-lateral motion is as- 
sumed. Thus, If B were to the left of the fault, 
then the angle a would decrease in size with 
time. If B were to the right of the fault, the 
angle would Increase. Note that the situation 
is reversed if left-lateral motion is assumed. 

Once the location of the active faults were 
Incorporated into the model, the data was read- 
justed several times. The first readjustment 
revealed an Inadequacy of the model in that It 
could not accommodate large discontinuities In 
motion as a function of time such as those which 
occurred along the Imperial fault as a result of 
the 1940 El Centro earthquake. To continue with 
the test all pre-1940 data was removed from the 
data set except for a 1935 astronomic azimuth 
near station YOTE (see figure 1). This azimuth 
was retained for better orientation control over 
time. In retaining the azimuth it Is assumed 
that the distance between the location of this 
observation and that of the earthquake Is suf- 
ficiently large that the orientation of the 
observed line did not change dlscontlnuously at 
the time of the earthquake. Future mathematical 
models need a feature to accommodate the rela- 
tively instantaneous shifts in position associ- 
ated with major earthquakes. 

With the pre-1940 data removed, the remaining 
data revealed a velocity pattern corresponding 
to a general shrinking of the network. Indicat- 


ing a problem with scale. The 1959 highway 
traverse, the 1967 observations, and the two 
Transcontinental Traverse projects 1974-76 all 
have sufficient electronic distance measurements 
to render good scale control for these epochs. 
However, closer examination of the 27 observed 
distances of 1959 indicated that they were too 
long by an estimated 13 parts per million. 

These observations were accordingly rescaled 
for the purpose of this test. Further investi- 
gation of these distances is being undertaken. 

An adjustment with the rescaled distances 
revealed a velocity pattern corresponding to a 
rotation about the fixed station, indicating an 
orientation problem. The data includes 13 as- 
tronomic azimuth observations, one observed in 
1935, one In 1967, and the remaining 11 as part 
of the two Transcontinental Traverse projects 
1974-76. These azimuths are not suspected to 
contain any serious non-random error. Instead 
the results Indicate a case of modeling obser- 
vational errors as movement. The model is a 
second degree polynomial in time and the above 
azimuths essentially represent three epochs, 
i.e., three points on the graph of network ori- 
entation versus time, the minimum required to 
determine a second degree polynomial. Hence 
all error in these three epochs of orientation 
is absorbed into the model. For the first two 
epochs, 1935 and 1967, the orientation is de- 
termined by a single observation. It is not 
unreasonable for an astronomic azimuth to be 
in error by one arc second which corresponds to 
0.485 meters at distance of 100 km from the 
fixed station. Further evidence of this effect 
was revealed by the high correlation coeffi- 
cients between the estimated parameters which 
are linear in time and the corresponding param- 
eters which involve the second power of time. 
Consequently, the data does not allow for the 
solution of an acceleration term. One way to 
overcome this weakness in the data would be 
to enlarge the network so as to Include addi- 
tional orientation control from nearby projects. 
There is a limit, however, to the effectiveness 
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Fig. 2. The position of station B relative to 
the fault can be established by the direction 
of change in the angle a with time. 
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of this technique. The more distant the ob- 
servations are from the location of interest, 
the smaller is their relative information 
content. 

With the elimination of all terms involving 
the second power of time, a solution involving 
16 time parameters was obtained. In this solu- 
tion the polynomial in the ith area is of the 
form 


^i.j = ‘’l,j,l +‘’i,j,2 




where l<i<3, l<j<2 and ($,X) are the coordinates 
assigned to the station ORIENT. The constraints 
^1 1 were Imposed. This corre- 

sponds to tAe’ assumption that station ORIENT in 
subregion 1 did not move with time. The esti- 
mated values for several parameters in this 
solution were below the estimated values of 
their sta.ndard errors. Hence, additional con- 
straints need to be Imposed to compensate for 
the inadequacies of the data. Some experimental 
adjustments were performed constraining differ- 
ent combinations of weakly determined parameters 
to specific values. Figure 3 Illustrates the 
velocity vectors relative to station ORIENT 
obtained in one of these experiments. Here 
the five constraints bj 2 2 =b 2 i i “b 3 jlj 3 = 
b3^2 1 “ ^3^2,3 “ 0 were imposed* in addition to 


fixing station ORIENT. The heavy wavy lines 
in figure 3 correspond the subregion boundaries 
input to the solution and dividing the region 
into three subregions. The error ellipses in 
figure 3 indicate the 95% confidence limits for 
the velocity vectors. Note that error ellipses 
are relative to the origin and depend on the 
choice of constraints. 

Statistical analysis in the form of an F-test 
indicate that the 11 parameter solution of 
figure 3 is overconstrained relative to the 16 
parameter solution at the 0.01 significance 
level. A few more adjustments were attempted 
to find the optimum set of constraints utiliz- 
ing the statistical concept of fixing a para- 
meter whenever there is insufficient information 
to significantly estimate its value. Sometimes 
more realistic constraints can be derived from 
the physical theory itself. Both the 16 and 11 
parameter solutions result in nonsymmetrlc 
strain matrices for each subregion. Physically 
this corresponds to a rotation of the network 
with time. The average rotation of the network 
obtained from these solutions is of the order 
of 0.1 (10“®) radians/yr. Since the estimated 
standard error in astronomic azimuths is 5.3 
(10“®) radians [Strange and Pettey, 1977], this 
rotation is probably only noise in the 13 ob- 
served azimuths. If it is physically plausible 
that the overall network does not rotate with 



Fig. 3. The velocity vectors relative to station ORIENT as obtained 
from an 11 parameter solution of the 1941-76 data. The ellipses are 


the 95% confidence limits, and the heavy wavy lines represent the 


subregional boundaries supplied in the solution. 
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time, then a more realistic type of constraint 
is to restrict some or all of the strain 
matrices to a symmetric form. Experiments with 
this type of constraint could not be performed 
in time to include them in this report. 

Evaluation of the Model 

The results of the various solutions trans- 
late into north-south contraction in the Impe- 
rial Valley in addition to right-lateral strike 
slip motion along the faults. The interpreta- 
tion of north-south contraction was at first 
questioned as a possible result of poor scale 
control in the data. However, the same inter- 
pretation was obtained [Savage et al., 1978] 
from trilateration networks observed by the U.S 
Geological Survey between 1972 and 1978. The 


data used by Savage et al. and the data of this 
study have no observations in common. Possible 
corroboration of the north-south contraction 
hypothesis is also provided by the observed 
subsidence in the area to the southeast of the 
Salton Sea (vicinity of station CALIPATRIA in 
figure 1) relative to the area just north of 
the Mexico-California boundary [Reese, 1977] . 

This subsidence is based on three epochs of 
leveling data spanning the period 1972-76. 

To further check the accuracy of the model, 
the 16 parameter solution is compared in table 
1 to the results obtained by Miller et al. 

[1970] and Savage et al. [1978]. In general the 
velocity vectors of this study have a more north- 
south trend than those of the other two studies. 
However, the difference is not statistically 
significant. Since a symmetric strain matrix was 


TABLE 1. Comparison of velocities relative to OLD BEACH deduced from 1941-67 
triangulation [Miller et al., 1970], 1972-78 trilateration [Savage 
et al., 1978], and the 1941-76 triangulation and trilateration of 
this study. 


Station Epoch yj (East) P 2 (North) 

(mm/yr) (mm/yr) 


BUTTE 

41-67 

16 



2 



72-78 

3 

+ 

9 

-2 

± 12 


41-76 

2 

± 

8 

-15 

± 5 

OROCOPIA 

41-67 

1 



1 



72-78 

6 

+ 

4 

-4 

± 4 


41-76 

4 

+ 

9 

-8 

± 8 

ALAMO 

41-67 

0 



3 



72-78 

1 

+ 

2 

-1 

± 3 


41-76 

0 

+ 

3 

5 

± 3 

SODA . 

41-67 

-16 



9 



72-78 

-3 

+ 

3 

13 

± 4 


41-76 

-10 

+ 

5 

18 

± 10 

KANE 

41-67 

-20 



19 



72-78 

-7 

+ 

5 

16 

± 4 


41-76 

-10 

+ 

7 

23 

± 8 

FISH 

41-67 

-15 



26 



72-78 

-12 

+ 

5 

24 

± 4 


41-76 

-11 

+ 

10 

33 

± 12 

DIXIE 

41-67 

-11 



33 



72-78 

-22 

+ 

9 

27 

± 4 


41-76 

-13 

± 

14 

35 

± 9 

OFFSET 225 

41-67 

-5 



43 



72-78 

-31 

+ 

14 

22 

± 10 


41-76 

-15 

+ 

18 

38 

± 9 

CARRIZO 

41-67 

-15 



33 



72-78 

-19 

+ 

8 

. 29 

± 5 


41-76 

-13 

+ 


36 

± 13 

OFFSET 229 

41-67 

-16 



35 



72-78 

-24 

+ 

14 

33 

± '5 


41-76 

-16 

+ 

19 

40 

± 12 


9J. 



assumed by Savage et al., the standard errors 
associated with their values are in general 
smaller than those of this study. 

In table 1 the numbers representing the solu- 
tion of this study would change by varying the 
number of parameters which are estimated. Some 
variations of the model that deserve investi- 
gation are the inclusion of the Sand Hill fault 
which runs from station GRAY northwesterly to 
station FRINK (figure 1) , and the inclusion of 
a fault line in the vicinity of station OFFSET 
227. In addition to varying the allowable coef- 
ficients of this model, it is desirable to re- 
fine the entire mathematical model so that the 
estimated parameters correspond closer to phys- 
ically observable quantities like the elasticity 
of the crust or the depth of faulting. On the 
other hand, even if the model were perfect the 
results of this study could differ from the 
results of Miller and Savage. The results of 
Miller are based only on the 1941 and 1967 net- 
works with the assumption that three stations 
were fixed in time. The results of Savage rep- 
resent different data over a significantly 
shorter time span and were obtained by assuming 
one station and one azimuth fixed in time. Fi- 
nally, recall that the model is unable to ex- 
tract acceleration information from this particu- 
lar data set . Thus it is impossible to check 
Thatcher's [1978] result that the average veloc- 
ity across the extent of the fault zone decreased 
from 82 ± 11 mm/yr for the 1941-54 period to 
23 ± 15 mm/yr for the 1954-67 period. 

Conclusion 

The pilot teat in the Imperial Valley demon- 
strated the advantages of fitting a model to 
the data. ' In particular, data from several 
sources can be assimilated. Observations of 
scale, orientation, and trlangulatlon which 
could not be used directly by the technique of 
comparing two sets of observations over the 
same quantities have been Included in a single 
data set. In the same manner the model will 
allow for the merger of classical geodetic ob- 
servations with data derived from radio inter- 
ferometry, creepmeters, Doppler, and satellite 
lasar-ranglng observations. However, before 
embarking on such an ambitious project, a model 
is sought which corresponds more closely to 
physical reality. This pilot test was a pre- 
liminary step in constructing such a model. It 
provides a departing point for future models and 
it reveals to some extent the Information con- 
tent of classical geodetic data. 
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